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We give a rigorous and mathematically well defined presentation of the Covariant and Gauge 
Invariant theory of scalar perturbations of a Friedmann-Lemaitre-Robertson- Walker universe for 
Fourth Order Gravity, where the matter is described by a perfect fluid with a barotropic equation 
of state. The general perturbations equations are applied to a simple background solution of ii" 
gravity. We obtain exact solutions of the perturbations equations for scales much bigger than the 
Hubble radius. These solutions have a number of interesting features. In particular, we find that for 
. all values of n there is always a growing mode for the density contrast, even if the universe undergoes 

' an accelerated expansion. Such a behaviour does not occur in standard General Relativity, where 

I as soon as Dark Energy dominates, the density contrast experiences an unrelenting decay. This 

■ peculiarity is sufficiently novel to warrant further investigation of fourth order gravity models. 

PACS numbers: 04.50.+h, 04.25.Nx 

I. INTRODUCTION 

Q^, In the last few years the idea of a geometrical origin for Dark Energy (DE) i.e. the connection between DE and a 

■ non-standard behavior of gravitation on cosmological scales has attracted a considerable amount of interest. 

(2jr)| Higher order gravity, and in particular fourth order gravity, has been widely studied in the case of the Friedmann- 
'— Lemaitre- Robertson- Walker (FLRW) metric using a number of different techniques (see for example [H, H, 3, B ) ■ 
_q ' Recently a general approach was developed to analyze the phase space of the fourth order cosmologies [3, llOt . 
^ , providing for the first time a way of obtaining exact solutions together with their stability and a general idea of the 
\^ ' qualitative behavior of these cosmological models. 

, The phase space analysis shows that for FLRW models there exist classes of fourth order theories which admit a 
transient decelerated expansion phase, followed by one with an accelerated expansion rate (see also [7| for a different 
approach). The first (Friedmann-like) phase provides a setting during which structure formation can take place, 
followed by a smooth transition to a DE-like era which drives the cosmological acceleration. However, in order to 
' determine if this is the case, we need to develop a theory of cosmological perturbations for higher order gravity. 
, The aim of this paper is to give a rigorous and mathematically well defined theory which describes the evolution 
' scalar perturbations of FLRW models in fourth order gravity, which can be used to investigate this issue in detail (see 
pH ll% [Tsl for other recent contributions to this area) . 
In order to achieve this goal a perturbation formalism needs to be chosen that is best suited for this task. One 
. , possible choice is the Bardeen metric based approach [l^ which guarantees the gauge invariance of the results. 
. However this approach has the drawback of introducing variables which only have a clear physical meaning in certain 
" " gauges [l^. Although this is not a big problem in the context of General Relativity (GR), this is not necessarily true 
in the case of higher order gravity and consequently could lead to a miss-interpretation of the results. 

In what follows we will use, instead, the covariant and gauge invariant approach developed for GR in [TtI. [isl. [Tol. 
[20l . [2H which has the advantage of using perturbation variables with a clear geometrical and physical interpretation. 
Furthermore, we use a specific recasting of the field equations that will make the development of the cosmological 
perturbation theory even more transparent. 

The outline of the paper is as follows. In Section II we give the general set up of the equations. In Section III we 
introduce the 1-1-3 covariant approach and apply it to the field equations presented in Section I. In Section IV we give 
the general evolution equations for the 1-1-3 quantities and linearize them around a FLRW spacetime. In Section V 
we introduce the perturbation variables and their propagation equations giving also their harmonic decomposition. 
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In Section VI we apply these equation to the case of i?"-gravity with a barotropic perfect fluid matter source and find 
exact solution in long wavelength limit. Finally in Section VII we present our conclusions. 

Unless otherwise specified, natural units {h ~ c = ks = SttG = 1) will be used throughout this paper, Latin indices 
run from to 3. The symbol V represents the usual covariant derivative and d corresponds to partial differentiation. 
We use the —,+,+,+ signature and the Riemann tensor is defined by 

R\cd = W\d^a - W\c4 + W\dW\, - WhcW%f , (1) 

where the are the Christoffel symbols (i.e. symmetric in the lower indices), defined by 

Wt^d = ^g"" (.9be,d + 9ed,b - gbd,e) ■ (2) 

The Ricci tensor is obtained by contracting the first and the third indices 

Rab — g'^'^Racbd ■ (3) 

Finally the Hilbert-Einstein action in the presence of matter is given by 

(4) 



dx 



II. GENERAL EQUATIONS FOR FOURTH ORDER GRAVITY. 

The most general action for a fourth order theory of gravity is given by 



~g[K + coR + ciR^ + C2Rt,,R^'' +CJ\ , 



(5) 



where we have used the Gauss Bonnet theorem [23 ] and Cm represents the matter contribution. In situations where 
the metric has a high degree of symmetry, this action can be further simplified. In particular, in the homogeneous 
and isotropic case the ^ can be expressed as 



A = j d^x^ifiRj + Cn] 



Varying the action with respect to the metric gives the generalization of the Einstein equations: 



f'Gab — f ( Rab — 2 gabR 
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gab {R - Rf) + ^b^af - SafcV.V^/' 



(6) 



(7) 



where / = f{R), f - 



andT-. 



2 Siy^Cra) 



represents the stress energy tensor of standard matter. 



dr ' " V-5 Sg^il' 

These equations reduce to the standard Einstein field equations when f{R) = R. It is crucial for our purposes to be 
able to write ([7]) in the form 



Gab = T 



ab 



rr>H rptot 

-'^ab ^ah ' 



where T, 



ab 



^ab 
f 



and 



rpR. 

^ab 



1 

7 



-(Jab (R - Rf) + ^b^af - gab^c^'f 



(8) 



(9) 



represent two effective "fluids": the curvature "fluid" (associated with T^) and the effective matter "fluid" (associated 
with T^). This step is important because it allows us to treat fourth order gravity as standard Einstein gravity in 
the presence of two "effective" fluids. This means that once the effective thermodynamics of these fluids has been 
studied, we can apply the covariant gauge invariant approach in the standard way. 

The conservation properties of these effective fluids are given by the Bianchi identities T^J^* . When applied to the 
total stress energy tensor, these identities reveal that if standard matter is conserved, the total fluid is also conserved 
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even though the curvature fluid may in general possess off-diagonal terms 123, l2J| . In other words, no matter how 
complicated the effective stress energy tensor Tj^* is, it will always be divergence free if T™''' = 0. When applied to 
the single effective tensors, the Bianchi identities read 

irn;b r/f 

C R-' , (10) 

(11) 

with the last expression being a consequence of total energy-momentum conservation. It follows that the individual 
effective fluids are not conserved but exchange energy and momentum. 

It is worth noting here that even if the energy-momentum tensor associated with the effective matter source is not 
conserved, standard matter still follows the usual conservation equations = 0. It is also important to stress that 

the fluids with and defined above are effective and consequently can admit features that one would normally 
consider un-physical for a standard matter field. This means that all the thermodynamical quantities associated with 
the curvature defined below should be considered effective and not bounded by the usual constraints associated with 
matter fields. It is important to understand that this does not compromise any of the thermodynamical properties of 
standard matter represented by the Lagrangian Cm ■ 



III. COVARIANT DECOMPOSITION OF HIGHER ORDER GRAVITY 



In this section we will describe the general covariant decomposition fourth order gravity. This procedure will take 
place in two steps. The first one is to develop the kinematics of the spacetime and the second one is to study the 
thermodynamics of the effective fiuids defined in the previous section. 



A. Preliminaries 



The starting point for any analysis using the covariant approach is the choice of a suitable frame i.e. the 4- velocity 
Ua of an observer in spacetime. 

This choice changes the structure of the equations and can simplify the calculations in the same way a choice of 
coordinates makes life easier in classical mechanics. Even if the covariance of the theory guarantees that all velocity 
fields are equivalent, a number of natural choices for Ua exist; they are the energy frame , which is defined to be a 
timelike eigenvector of the stress energy tensor Tab, the particle frame that is derived from the particle fiux vector 
Na and the entropy frame wf defined by the entropy flux vector Sa- These frames have the advantage of inducing 
important simplifications to the equations (e.g. in the energy frame the total energy fiux is always zero [T6|l. 

The vectors u^, , always exist in the case of a perfect fluid and coincide, defining a unique hydrodynamical 
4- velocity for the fiow [l^. In the case of more than one perfect fluid, one can in principle define these three frames 
for each component as well as for the total fiuid and then choose the most convenient frame to work in. 

In our specific case, equation ([5]) allows us to define two "effective" fluids, but their structure does not necessarily 
make any of the three frames defined above a suitable choice. This because both the curvature fluid and the effective 
matter do not necessarily satisfy the Weak Energy Condition (WEC) TabVV'' > 0, VVa < (Hi. This relation, 
which for homogeneous and isotropic spacetimes corresponds to the requirement that the energy density is positive, 
is the key hypothesis which allows the timelike vectors Mf',M^,Mf to exist and is, in general, a very reasonable 
assumption. In our case, however, the violation of this condition means that none of the energy entropy or particle 
frames are, in general suitable choices of frame. 

However, an alternative frame choice follows from the fact that whatever the behavior of the effective fluid is, 
standard matter is still thermodynamically well defined and consequently the stress energy tensor Tab satisfies the 
standard energy conditions. It follows that a natural choice of frame the one of those associated with standard matter 
{ua = w™), assumed to be a barotropic perfect fluid with equation of state p — w^. (33 |. This choice is also motivated 
by the fact that the real observers are attached to galaxies and these galaxies follow the standard matter geodesies. 
Consequently this frame choice is the one which can be best motivated from a physical point of view. 
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B. Kinematics 

Once the frame has been chosen the derivation of the kinematical quantities can be obtained in a standard way [2^ . 
The derivative along the matter fluid flow hues is defined hy X = UaV^X . The projection tensor into the tangent 
3-spaces orthogonal to the flow vector is: 

hah = gab + UaUb /l^fe/l^ = /l°c , habU^ = , (12) 

and the tensor VhUa can be expanded as 

V&Ua = VfeWa - ttaUb , ^ bUa = + CFab + ^q6 , (13) 

where Vq — h^a^b is the spatially totally projected covariant derivative operator orthogonal to w". This relation 
allows us to define the key kinematic quantities of the cosmological model: the expansion 0, the shear aab, the 
vorticity ujab and the acceleration Ua — Ua- 

In the following, angle brackets applied to a vector denote the projection of this vector on the tangent 3-spaces 

V(a) = ha'Vb . (14) 

Instead when applied to a tensor they denote the projected, symmetric and trace free part of this object 

W^ab) = [ha'h)" - ^h-'^hab] ■ (15) 

Finally the spatial curl of a variable is 

{xy^ = if^^" VcX''>d (16) 

where tabc — u'^Vabcd is the spatial volume. 

The general propagation equations for these kinematic variables, for any spacetime corresponds to the so called 
1+3 covariant equations ^] and are gven in Appendix [Cl 



C. Effective total energy-momentum tensors 

The choice of the frame also allows us to obtain an irreducible decomposition of the stress energy momentum tensor. 
In a general frame and for a general tensor Tab one obtains: 

Tab = llUaUb + phab + ^qi^a^b) + T^ab , (17) 

where /x andp are the energy density and isotropic pressure, qa is the energy fiux (qa — q{a)) a-nd TTab is the anisotropic 
pressure {-Kab = Tr{ab))- 

This decomposition can be applied to our effective energy momentum tensors. Relative to u™ we obtain 

^tot ^ T^'°tw''u^ = = ^T'^Th"'' - p"'+p" , (18) 

O 

„tot T^tot i^h ~m I ^tot T^tot uc i-,d 2Lm \ „ R ( ac\\ 

la = --^bc "-a" = Qa + Qa : ^ab = ^ cd "'Ka'^by = ^ab + ^ab ^ (19) 

with 

m m m m 

~m P ^ -m ^a ~ m ' ab /'0^^ 

Since we assume that standard matter is a perfect fluid, q™ and ttJ^ are zero, so that the last two quantities above 
also vanish. 
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The effective thermodynamical quantities for the curvature "fluid" are 

"1 
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(Rf - /) - ef'R + f'V'R + f" iibVR 



\{f- Rf) + f'R + V"'R' + l^fR If'^'R- 
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f"'R\IaR+ f'VaR- - f'^aR 
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f"V^a^t)R+f"'^iaR^b)R + <yabR . 

The twice contracted Bianchi Identities lead to evolution equations for /i™, /i^, q^: 



/'2 ' 



R, , ,,nif"^aR, 



(21) 

(22) 
(23) 
(24) 



(25) 
(26) 



qfa) + VaP^ + V^f, = - I e - - (^^ + P^) ~ " ^6 + ^"^^^^T^^eq : consZ , (27) 

and a relation connecting the acceleration iia to /i™ and p™ follows from momentum conservation of standard matter: 

V'^p" = -(^™ +p")mV (28) 

Note that, as we have seen in the previous section the curvature fluid and the effective matter exchange energy and 
momentum. The decomposed interaction terms in Equations (j26p and (|25p are given by ^™ f J''^^ Z^™ ^ ■ 



IV. PROPAGATION AND CONSTRAINT EQUATIONS 
A. Nonlinear equations 

We are now ready to write the full nonlinear 1+3 equations for higher order gravity. Substituting the quantities 
given above into the equations given in the Appendix A, we obtain the following results: 



Expansion propagation (generalized Raychaudhuri equation): 

e + + aab'j"''' - 2u:aUj"- - V"-Ua + Uail" + \{pJ^ + 3p™) = -^(a*^ + 3p^) , 

Vorticity propagation: 

LO(a) + |0^a + ^CUrlUa - (7ab^^ = , 

Shear propagation: 
Gravito-electric propagation: 

E(^ab) + QEab - CUrliJab + ^(m™ +P")crafc " 2?i'=£crf(a-f^fc) - 3crc{a-Bh)'= + aj'^e^diaEb)' 

Gravito-magnetic propagation: 

H(ab) + 0-f^afc + CMvlEab — ^i^^Hb)'^ + ^'^ ^ cd{aE[ b)'^ + '^u'^^ cd{aEb)'^ 



(29) 
(30) 

(31) 
(32) 
(33) 
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Vorticity constraint: 



Shear constraint: 



Gravito-magnetic constraint: 



V^iUa - il^'iOa = , (34) 
VVab - curlw, - fV^e + 2[LU,u]a = , (35) 

curl (Tafc + V(aWt) - Hab + '2U{a^b) = , (36) 



Gravito-electric divergence: 

^'Eab - i V,m" - + SHabCo" = ia^qf - | [c^, q\ ~ i VVf, + i v,^« - iOgf , (37) 

Gravito-magnetic divergence: 

\/'Hab - ir + Pn^a + [^7, E]a - S^afcC^" = -^CUrlgf + (^i^ + p'')uJa " ^ ^ 7r«], - ^Tr^V^ . (38) 

The standard GR equations are obtained by setting /(i?) — R which corresponds to setting aU the right hand sides 
of these equations to zero. Together with Eqs. these equations govern the dynamics of the matter and 

gravitational fields in fourth order gravity. As we will see the new source terms in the propagation and constraint 
equations will modify the evolution of the perturbations in a non-trivial way. 

B. Linearized equations 

In the previous section we derived the exact nonlinear equations that govern the exact gravitational dynamics of 
fourth order gravity relative to observers comoving with standard matter. These equations are fully covariant and hold 
for any spacetime. Consequently, we can linearize these equation around any chosen background, avoiding the need for 
choosing coordinates and dealing directly with physically well defined quantities, rather than metric components [ssj . 
These features, which are desirable in the GR case, become essential for the correct understanding of the evolution 
of perturbations in fourth order gravity as well as in other kinds of alternative gravity theories [l^l • 

In what follows we will choose a Friedamnn-Lemaitre-Robertson- Walker (FLRW) metric as our background. We 
make this choice for a number of different reasons. First of all the possibility of writing a general fourth order 
lagrangian as a simple function of the Ricci scalar is surely possible for this metric. Secondly, because most of the 
work in GR perturbation theory has been performed for this background it makes a comparison of behavior of GR 
and fourth order gravity more straightforward. 

The Friedmann background is characterized by the vanishing of all inhomogeneous and anisotropic quantities , tt^j, 
and defines the order of the quantities appearing in the 1-1-3 equations and the linearization procedure. In particular, 
the quantities that are zero in the background are considered first-order of in the linearization scheme. In addition, 
the Stuart & Walker lemma ensures that since these quantities vanish in the background, they arc automatically 
gauge-invariant . 28.] . 

The cosmological equations for the background read: 

= 3/i™ + 3^^ - I , (39) 
e + -f i(/i'" + 3p") + i(/i^ H- 3p^) = , (40) 

A™ + e (Ai" + p") = , (41) 

where ji^ and are the zero order energy density and pressure of the curvature fluid, 7? is the 3-Ricci scalar and 
R = QK/S^ with the spatial curvature index K = 0,±1. 
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Linearization of the exact propagation and constraint equations about this background then leads to the system: 



e + ie^ - vM, + + 3p™) = + 3p^) , (42) 

LUa + 2HLUa + |curl = , (43) 

&ab + '^HUab + Eab - V(aA) = -Qa , (44) 

Eab + SHEab - cnrl Hab + i(/i" + P")a,fa 

- -Uf^^'+p'^Wb - ^irf^b) - l^iaQb) ~ l^^ab , (45) 

Hab + iHHab + curl Eab - jcurl^ffc , (46) 

VVab - curl LOa " | VaO - -gf , (47) 

CUrldaf, + V(at^h) - i?ab = , (48) 

y^Eab - iVaA" = -^V^TT^, + iVaM"" ^ ^©9^ , (49) 

V^Hab - (A" + P^UJa = -^CUrlgf + (m"" + P^'Oc^a , (50) 

V'^yja - , (51) 
together with the linearized conservation equations 

/i" = -e(/"+p'") , (52) 

VV" = -(Ai^+p™)!*" , (53) 

M« + V"qf - e (m^ + P^) + ^ , (54) 

9fa) + V.p« + V^f, = - I e - (/.^ + p«) + A^™^^^^ , (55) 



obtained from ((25)) - ([28l) . Note that at first order the equation of the vorticity ((5T|) is homogeneous i.e. the evolution of 
the vorticity is decoupled. This will be important in the next section when we will derive the perturbations equations. 
These equations provide the basis for a covariant and gauge-invariant description of perturbations of f{R) theories of 
gravity. 

V. DYNAMICS OF SCALAR PERTURBATIONS 
A. Perturbation Equations 

We are now ready to analyze the evolution of the density perturbations on a FLRW background. The quantities 
appearing in the linearized equations given in the previous section can be decomposed in scalar vector and tensor 
components, i.e. 

= Va + K = ?y"'"Vf,Fc + V"?) , where = , r;'^''=VbT4 = , (56) 

and 

Wab = Wab + Wab + T^„*fc = Wab + V aWb + V aV bW* , (57) 

where 

V"W'a6 = 0, iW)ab = 0, iW*)ab = 0, (58) 

and both of these decompositions are unique. Note that here we define scalars, vectors or tensors as quantities 
that transform like scalars, solenoidal vectors or symmetric tensors, or are obtained from them using the hab or Vq 
operators [181] . 

In linear regime and in homogeneous and isotropic backgrounds these different components do not interact with 
each other. In the following we will focus only on the evolution of scalar perturbations because they are directly 
related with density fluctuations. This can be done simply discarding the non scalar quantities in the equation above 
i.e. setting 



Va = VaV, Wab^V(^a^b)W. 



(59) 
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The identities in Appendix [Cl the vorticity constraint equation (I34p and the gravito-magnetic constraint equation 
([55)1 then show that 



CurlK = = CnrXWab , V''VKah - ^V^{VaW) , Ua^O^Hab, 



(60) 



as in standard General Relativity. 

In order to derive the equations governing density perturbations in the general case, we define the density and 
expansion gradients 



a ,,m "'^ ' 



Ca — SWaR 



and the (dimensionless) gradients describing inhomogeneity in the Ricci scalar: 



(61) 



(62) 



Another important quantity in the treatment of the evolution of the density perturbations is the Newtonian potential 
(defined through the divergence of the electric part of the Weyl tensor ([57]) [IS]). 
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29/' + 3i?/" 

7^a 



(63) 



where represents the total energy density fiuctuation. 

Using equations equations the identities in Appendix [Cl a ssuming matter to be a barotropic 

perfect fiuid with barotropic factor w = p^Z/x™ and that the vorticity is zero [36|, we obtain the following system of 
evolution equations for the above variables: 



{l+w)Za 



Rf" 29 
~1' 
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3(w - 1) (3w + 2) ^ 2^92 + iw{ii" + 3p^) 
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9 + 2R^— 1 - RZa 



w 



3 /" w + 1 /" 3(w + 1) 



i? 9 - 3i? 



if 



52 

+V27^a , 

together with the constraint 



i/: + ^7?2 + 9^7?-^92 + V + 3p-)+i?:^-U + i(.. + l)ii 
3/" /' /' 9 ^3V^ ^ j„ 6/' 2^ V 
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(64) 

(65) 
(66) 

(67) 



52 



+ u© 



2i?/" 



2i?9^ - ^ ( / - 2/i + 2i?9/" + 2 



7^a 



29/" 



2f" ~, 



/' 



(68) 



9 



The propagation equation for the variable C is 

36/"7^a 36/'!?™ 



Ca = 



S'2 (29/' + 37?/") 5-2 (26/' + 3Rf' 



6/' 



16^6 4/'e2 - 12/'/i^ 



" \3{uj + l) 2e/' + 3i?/" J 2e/' + 3i?/ 



5*2 ( 26/' + 3Rf" 
120 f" 



Ca 



29/' + 3Rf" /' 



2^7 )3?a 



25*2 fe/" - 37?/(3)) 127?e/'/(3) - 2/" (3/ - 2 (62 - 3^*^) /' + 6Ref" 



3/' 



3(c^+l) " /' 



^29/' + 3i?/" ) /' 
2^2 [e/"_3i^/(3)'j 

3/' " 



(69) 



this equation, which is redundant, will be used in Section IVTl to substitute (|66p because of its specific form in the long 
wavelength limit [2^. 



B. Scalar Variables 



The variables we have defined above describe the general evolution of the density perturbations and the other 
scalars on a FLRW background. The phenomenon of the clustering of matter is traditionally described, however, 
considering only the scalar part of these variables. This can be easily done using the local decomposition [I^l 



S\7°'Xa = Xab = -^habX + S^j, + X[ab] 



where 



^ab = X{ab) - -^habX . 



(70) 



so that the operator Vq applied to the (|6ip and ((62)) extracts the scalar part of the perturbation variables. In this 
way we can define the scalar quantities 



.AT 



(71) 



which will characterize the evolution of the spherically symmetric part of the gradients (|61ll62p . The evolution 
equations for the first four of these variables are 



A™ = w9A,„-(l + w)Z, 



S{w - l)(3w; + 2) fi 2wQ^ + 3w(/i^ + 3p^) 
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(75) 
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together with the constraint 
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In standard GR, only the first two equations and the last one are present and the density perturbations are governed 
by a second-order equation for A™ whose independent solutions are adiabatic growing and decaying modes. The 
presence of fourth order corrections introduces important changes to this picture. In fact, in this case the evolution of 
the density perturbations is described by a closed fourth order differential equation which can be obtained form the 
above first order equations. This follows clearly from our two effective fluids interpretation. 



C. Harmonic analysis 



The system ([72l) -(|75 |) is a system of four partial differential equations which is far too complicated to be solved 
directly. For this reason, following a standard procedure we perform an harmonic decomposition. This allows one to 
reduce equations (|72[) - (|75p to ordinary differential equations which are somewhat easier to solve. 

In the covariant approach the harmonic decomposition is performed using the trace-free symmetric tensor eigen- 
functions of the spatial the Laplace-Beltrami operator defined by [l6| : 

v^g = -^g , (78) 

where k = 27rS'/A is the wavenumber and g = 0. Using these harmonics we can expand every first order quantity in 
the equations above [st! ]. 

X(i,x) = ^xW(t)gW(x) (79) 
where ^ stands for both a summation over a discrete index or an integration over a continuous one. 
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Developing in terms of Q, ^^-^8$ reduce to 
GAl^J - (1 + , 

3{w - l){3w + 2) ^1 2we'^ + 3wi^i^ + 3p^) w k 



3ff(fc) 



1 _ r_^_ 1 i^r _ 

2 2/'/' /'/' 





6{w 


+ 1) 


W+lS'2 






7^«, 











w + l 



-R 



e + '^Ryjrj ^^''^ - ^^^^^ 



3 /" w + l j" 3(iy + 1) 



i? ( e - 3i?^ 



+K 



36/'A('=) 



Tlik) 
12/" 



(80) 

(81) 
(82) 

m 

(83) 



52 (29/' + 3i?/") 52 (29/' + 3i?/") 5*2 (29/' + 37?/" 



6/' 



S'2 (29/' + 37?/' 



4cj9 4/'92 - 12/'/i^ 



129/'' 



29/' + 3i?/" 



+ 1 29/' + 3i?/" 
127?9/7(3) - 2/" (3/ - 2 (92 - 3^^^) /' + 6i?9/" 



29/' + 3i?/" /' 



k 



4u;S^e . 2^2^. ^^^^ 2^2 (0^. „ 3^^(3) 



3{uj + l) 



/' 



3/' 



(84) 



= 



52 



^e + ^zW-2-^AW + 



2i?9^ _ £_ (^/ _ 2/i + 2i?9/"J - 2^ 



/' /' 



/' ^2 



/' 



(85) 



Finally it is useful to write equations (|80ll83p as a pair of second order equation. In this way the GR limit is more 
transparent when the written in this form. They are: 



A(fc) - 


f 2\^^7?/" 

r'3j®+ /' 


A(fe) - 













^ 2ujRQf" (3tj2 _ ^ 



/' 



/' 



2(^ + 1) 



/" 



k 



K 



f'ii^''^ + (9/" + 2i?/(3)) TZ^f^) 

(/(3)i?2 + (p^ + ^^)/' + imf + i?/") 



f 

2 



7^(.) _ i-^ + mr ^jk) 



1 



6 



;{ZW - l)/i 



aw - i^-mf\ ik^ 

w + l 



In the GR limit we have / = i?, so the above equations reduce to 

K^k) _ _ I) eAW - [c,fc2 - (i + - a« - 

7^ = (3w - 1) a^ac^^) . 



(86) 



(87) 



(88) 
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The second of these equations is just the spatial Laplacian of the trace of the Einstein Field equations i? = 3p — /i. 



VI. EXAMPLE: i?"-GRAVITY 



Let us now apply the equations derived in the above sections to the simplest example of fourth order theory of 
gravity: _R"-gravity. In this theory f{R) — X-R" ^nd the action reads 



A = 



where x a the coupling constant with suitable dimensions and x 
If i? ^ the field equations for this theory read 



M\ 



1 for n ~ I. 



(89) 



Gab = X~ 



-_j_ab_ 



T 



where 



^ab — 



^ab — 



rpM 

-l_!_ab_ 



R 



+ (n - 2) 



i?2 



{gcaQcLb — gcdQab) 



(90) 

(91) 
(92) 



The FLRW dynamics of this model has been investigated via a complete phase space analysis in [s*]. This analisys 
shows that for specific intervals of the parameter n there is a set of initial conditions with non zero measure for which 
the cosmic histories include a transient decelerated phase which evolves towards an accelerated expansion one. This 
first phase was argued to be suitable for the structure formation to take place. 

In what follows we will analyze the evolution of the scalar perturbations during this phase in the long wavelength 
limit. In this approximation the wavenumber k is considered to be so small that the wavelength A = 2TTS/k associated 
with it is much larger than the Hubble radius. Equation (l78l) then implies that all the Laplacians can be neglected 
and the spatial dependence of the perturbation variables can be factored out. It is also well known [2^ that in this 
limit and in spatially flat {K = 0) backgrounds the reduces to C = i.e. the variable C is conserved so that the 
number of perturbations equations can be reduced to three. 

Let us now set the background to be the transient solution 



S'oi='(i+™) , 



, /I = fiot 



-2n 



(93) 



of ^. The expansion, the Ricci scalar, the curvature fluid pressure, the curvature fluid energy density and the effective 
matter energy density take the form: 



e = 
R = 



2n 



An[4n-3{uj + 1)] 
3i2(w + l)2 ' 

2(n-l)[2n(3w + 5)-3(cj+l)] 
3<2(t^ + l)2 ' 

2(n - 1) [n (6u;2 + 8w - 2) - 3cj(w + 1)] 



nx 



3t2(tJ+ 1)2 
n(4n-3(w + l)) 

<2(t^+l)2 



4n2 - 2(n - l)[2n(3w + 5) - 3(w + 1)] 



3(a 



l)2t2 



(94) 
(95) 
(96) 
(97) 
(98) 
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Substituting in the equations given above and passing to the long wavelength limit we obtain 
2n 6(n — l)n 



A 



w + 1 n + 3(n — l)w — 3 



\n 3(W+1)-^ ^ 1 4„ 

t 4a2[n + 3(n-l)w-3] ° 



3(„ _ + i)^[n{Qw + 8) - 15(^; + 1)] ^ ^ 9(n - l)(u; + 1)^^^ , 

4[n + 3(n-l)w-3][4n-3(?z; + l)] 4[n + 3(n - l)^ - 3] [4n - 3(w + 1)] ' ^ ' 

^^^^ 8,..,(4n-3(^ + l)) A^^ (100) 
3(w + Ij'^ t-^ 

■ 2n(4n-3w-3) r7n 4„ . ^/ 3n(n-l) n „ \ 

(u) + l)(n + 3(n-l)w-3) a2 \n + 3in - l)w - 3 w + 1 )t 

(Mn-2)(n I) _ 3n^(9"-26) + 57 8n^(n 2) _ \ U 

V " + 3(n-l)w-3 9(w + l)(n-l) 9(w + l)2(n-l) / V y 



16n(4n + 3(n - l)it; - 3)(4n - 3(w + 1)) {{9w{w + 1) + 8)n2 - (3w(9w + 8) + 13)n + 2,{w + l)(6u; + 1)) 

27(n - l)(w + l)4(n + 3(7i - l)w - 3) 



where Co is the conserved value for the quantity C . The evolution of density perturbations can then be decoupled 
via the third order equation 

(n-l)A„-(n-l) (±^-^\^+V,{n,w)^+V2{n,w)^+V^{n,w)C^t-^^-^=0 (102) 



where 



2(-9(2(n-l)n + lV2_[,6n(n(4n-7) + l)w + 18cj + n(4n(8n-19) + 33) + 9) , ^ 

= 9(i:tt? 

^^^^^ ^ 2((2n - 1)lo - l)(4n - 3(c^ + 1))(3(l^ + 1) + n(-9c^ + n{&LO + 8) - 13)) 

9(cj + 1)3 



7i(21w - 6n(cj + 2) + 31) - 18(w + 1) 
This equation admits the general solution 



A„ = Kit—^^ + K2t°'+ + /sTai"- - ^44^ ^iz;^ , (i06) 

at 



where 



1 nuj J{n - 1) (4(3tJ + 8)2^3 - 4(3w(18w + 55) + 152)^2 + 3(a; + l)(87w + 139)n - 81(t^ + 1)) 
q;± = - - H ± 



2 



2 uj + l 6(n-l)(a; + l) 

(107) 

9(tj + l)3(18(tj + 1) + n(-21cj + Qn{uj + 2) - 31)) 

~ 8(n(6w + 4) - 9(u; + 1)) (6(w + 2)^3 - (9w + 19)n2 - 3(cj + l)(3u; + l)n + 9(w + 1)2) ' *• ' 

Let us now focus on the case of dust {w — 0). The above solution becomes 

A„ = Xit-i+X2i"+'™="+i^3t"-'™=°-if4%i'~^ , (109) 

where 



I 1 , V(n-l)(n(32n(8n-19) + 417)-81) 
a±|_„ = ± ^^—^^ , (110) 

9(n(12n-31) + 18) 

A4k=o - 3^4^ _ ^^2n3 - 19n2 - 3n + 9) ' ^ ^ 

A graphical representation of the behavior of the exponent of the modes in (|109p as n changes is given in Figure [T] 
This solution has many interesting features. For 0.33 < n < 0.71 and 1 < n < 1.32 38] the modes i"±l™=o become 
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oscillatory. However since the real part of the exponents a±\uj=o is always negative the oscillation are damped and 
bound to become subdominant at late times. The appearence of this kind of modes is not associated with any 
peculiar behavior of the thermodynamic quantities in the background i.e. none of the energy condition are violated 
for the values of n which are associated with the oscillations. The nature of these oscillations is then an higher order 
phenomenon. Here we will not undertake a detailed investigation of the origin of these modes, such a study will be 
left for a future work. Also, for most of the values of n the perturbations grow faster in i?"-gravity than in GR. In 
fact only for 1.32 < n < 1.43 all the modes grow with a rate slower than i^/'^. 

Probably the most striking feature of the solutions p09p and (|106p is that the long wavelength perturbations grow 
for every value of n, even if the universe is in a state of accelerated expansion (see Figure[T]). This is somehow expected 
from the fact that in [H the fixed point representing our background is unstable for every value of the parameters. 
However, the consequence of this feature is quite impressive because it implies that in i?" gravity large scale structures 
can in principle also be formed in accelerating backgrounds. This is not possible in General Relativity, where it is 
well known that as soon as the deceleration parameter becomes positive the modes of the A solutions (or density 
contrast) are both decreasing. The suppression of perturbations due to the presence of classical forms of Dark Energy 
(DE) is one of the most important sources of constraints on the nature of DE itself. Our example shows that if one 
considers DE as a manifestation of the non-Einsteinian nature of the gravitational interaction on large scales, there 
is the possibility to have an accelerated expanding background that is compatible with the growth of structures. Of 
course, in order to better understand this effect, one should also analyze the evolution of perturbations on small scales. 
However this analysis is beyond the scope of this paper and it is left to left to a future, more detailed investigation. 

In the limit n 1 two of the modes of (|106p reproduce the two classical modes t^/'^ and typical of GR, but the 
other two diverge. At first glance this might be surprising but it does not represent a real pathology of the model. In 
fact equation (jl02p reduces to a first order differential equation when n = 1. Therefore in this case the two modes in 
the solution can be discarded and GR is recovered. 

From the system (l99t we can also obtain the solution for the other scalars: 

n = KrJ^-^ + KQtl^+ + Kji^^- - Kfi-^r^^ (112) 



= KqI^-^ + Kiot''+ + Kiif - - Ki2-^rT^(^-^ (113) 

where 



P± = a±-2, (114) 
7± = Q;± — 3 (115) 

and the constants K^, ..K12 are all functions of Ki, ..K4. These expression are rather complicated and will not be 
given here. It is interesting that these quantities have an oscillatory behavior for the same values of n for which A™ 
is oscillating. Also for these quantities the oscillating modes are always decreasing. 

Finally it is useful to derive and expression for the Newtonian potential <i>Ar given in (|63p which for our background 
takes the form 

Ana'^Kit^'^si^-^ 4„ (2ncj - {uj + l)a_) alK2t^i^'^'^- 

— — —TTT^ h 



3(t^ + l)2 3(cj + l)3 

4n {2nui - {uj + l)a+) alKst^i^^^^ 9{uj + if ~ 16n(2n + 3(n - l)w - 3)K4 
3(0; + 1)3 ^ 18(tJ + l)3 



Co ■ (116) 



As in the GR case, this potential has a constant mode and at least one monotonic mode. The presence of a constant 
mode it is important because it is consistent with the standard Sachs- Wolfe plateau in GR. In addition, the fourth 
order correction induce oscillations in two of the modes of (|116p as it is expected from (|63p and the form of the solution 
for Am- As for A^, these oscillatory modes always decay. Again, when n ^ 1 the Newtonian potential reduces to the 
one obtained in General Relativity. 
In the case of dust we have 



4n o ^, 4n , 4:71 r,^^ , in.f) I 4n , ^, , <lni« 



I- 1^(271 -3)K4 



(see Fugure[2]). This potential is weaker that the GR one at early times and became stronger at late times. The non 
constant modes are decaying only for —0.63 < n < 0.33, 0.71 < n < 0.86 and 1.32 < n < 1.36 and their exponent is 
grater than the one in GR for any n but 0.28 < n < 0.33. 
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FIG. 1: Plot of the real part of the exponents of each modes of the solution (|106|) against n in the case of dust {w — 0). 
The continuous and dashed line represent the modes t°'^ respectively (note how they coincide when a± are complex), the 

dashed-dot line represents the mode t 3(u,+i) ^nd the dot line the mode t 
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FIG. 2: Plot of the real part of exponents of the modes of the Newtonian Potential (??) against n in the case of dust (ui = 0). 
The continuous and dot-dashed line represent the modes i"3"+'3±l™=o respectively (note how they coincide when /3± are complex), 
the dotted line represents the mode t 3 . 

VII. CONCLUSION 

In this paper we have analyzed in a rigorous and mathematically well defined way the evolution of density pertur- 
bations of FLRW backgrounds in fourth order gravity. 

Our analysis has been based on two important steps. The first one follows from the fact that in homogeneous 
and isotropic spacetimes the field equations for a generic fourth order gravity theory can be rewritten in a form that 
resembles standard GR plus two effective fiuids. Then, using 1-1-3 covariant approach, it is possible to derive the 
general equations describing the evolution of the cosmological perturbations of these models in FLRW background. 
In this paper we have dealt only with the evolution of the scalar perturbations, the evolution of tensor and vector 
perturbations will be presented elsewhere [SOj, 31]. 

Providing that one has a clear picture in mind of the effective nature of the fluids involved, the approach above 
has the advantage of making the treatment of the perturbations physically clear and mathematically rigorous. In 
particular, it allows one to understand in a natural way that the equations governing the scalar perturbations in 
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fourth order gravity are of order four rather than two. 

Once the general perturbations equations were derived, we appHed them to the simple i?"-gravity model. In the 
long wavelength limit and using a background solution derived from an earlier dynamical systems analysis Q , we were 
able to find exact solutions to the perturbations equations. The results we obtain are particularly interesting. As 
expected, our background solution proved to be unstable under scalar perturbations - this solution always corresponds 
to a saddle point in the phase space of the i?" homogeneous and isotropic models . In addition, for specific intervals 
of values of the parameter n two of the four modes of the solution can become oscillatory and this might have 
consequences on the scalar perturbation spectrum. The connection between the spectrum of matter perturbations 
and the CMB power spectrum then offers an interesting independent way of testing these models on cosmological 
scales 

However the most striking property of the evolution of the density perturbations in this model is that growth 
is possible even if the background is accelerating. This means that unlike all the other models for dark energy, in 
R"' -gravity a decelerated phase is not necessarily required to form the large scale structure. Of course in order to fully 
support this claim a detailed analysis of the small-scale perturbations is necessary in order to understand better, if 
in this regime small scale structure formation is also preserved. However, the idea is very intriguing. A number of 
important constraints on the nature of Dark Energy come from the requirement that we have a matter dominated 
phase able to support the formation of large-scale structures. What we have discovered is that in i?"-gravity this is 
not necessarily the case. 

Also the Newtonian potential has a number of interesting features. First of all it contains a constant mode like in 
GR. This is an encouraging result as this is compatible with the Sachs- Wolfe plateau. On the other hand, the time 
evolution of the Newtonian potential is essentially different form GR. Such differences suggest that the dynamics of 
structure formation in i?"-gravity might be very different from what we find in GR and it is definitely worth a more 
detailed study. 

A natural question arises. How general are these results in terms of the form of the fourth order Lagrangian? Or in 
other words are also other, more popular, fourth order models able to give rise to the same effects? The question is not 
easy to answer. From the form of the decoupled equation for the matter density perturbations, we can conclude that 
features like the oscillating modes should be a common to all fourth order theories. On the other hand proving that 
density perturbation grow for all accelerating backgrounds is a much difficult matter. One interesting hint comes from 
the general dynamical system analysis given in ^9] in which backgrounds similar to the one we used for i?"-gravity are 
often unstable for most the values of the parameters of that theory. Finally, from the point of view of the dynamics 
of structure formation, it seems to us reasonable to say that if a model as "close" to GR as i?"-gravity has so many 
different features, fourth order gravity models with a Lagrangian very different from the one in the Hilbert-Einstein 
action will in general have a very different dynamics. However, as we have seen in i?" gravity, these differences do not 
necessarily imply a complete incompatibility with the data coming from the CMB and other observational constraints. 
However, much more work will have to be done before we can determine whether alternative gravity provides a viable 
alternative to our much cherished theory of General Relativity. 
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APPENDIX A: GENERAL PROPAGATION AND CONSTRAINT EQUATIONS 

For a general, imperfect energy-momentum tensor the propagation and constraint equations are: 
Expansion propagation (generalized Raychaudhuri equation): 

e + + CTafcff"'' - 2UaU'' - V'^ila + UaU" + 5 (M + 3p) 

Vorticity propagation: 
Shear propagation: 



ioia) + |BWa + ^CUrllia - CTah^^ = . (Al) 



^(ab) + fScTah + Eab - ^ {a^b) + CTciaf^b)" + ^^(0^6) - W(aMb) " l^T^ab = (A2) 

Gravito-electric propagation: 

+ ^'^{ab) + k^iaQb) + ^©TTab + ^d" (^a'^^^c + ^t{a'^b)d = (A3) 

Gravito-magnetic propagation: 

H{ab) + QHab + curl i?a6 - io^^^Hb)^ + Ld'^EcdiaHb)'^ + 2ifecd(aEb 



Gravito-magnetic constraint: 



ab) ~r 'TJIlab ^" cuii-C/afc - OO^/^a-^b) ^ ^cd{aJ^b) i^cd{a^b) 
= icurlTTab - lo^i^aSbyld (A4) 

Vorticity constraint: 

V°Wa - il'^Ua = . (A5) 

Shear constraint: 

VVab - CUrlcJa - fVaO + 2[W, U\a - = (A6) 
curl dab -I- V(at^6) " ^^ab + '^U(a^b) = . (AT) 

Gravito-electric divergence: 

^'Eab - iV.M™ - K i/]a + BiJ.bC^" - \alqb + |[^, g], + \V\ab - jVaA* + ieg^ = (A8) 
Gravito-magnetic divergence: 

^""Hab - {f^+p)0Ja + [<J,E]a - iEabUJ^ + ^CUrl^a + 5 k, Tl]a - ^TTafcCj'' = (A9) 

Here Ua = ^ea'"^wtic and the covariant tensor commutator is 

[W,Z]a=eacdW'eZ'''. 

APPENDIX B: COVARIANT FORMALISM VERSUS BARDEEN'S FORMALISM 

As we have seen the covariant approach is a very useful framework for studying perturbations in alternative theories 
of gravity. However, since most work on cosmological perturbations is usually done using the Bardeen approach [l5[, 
we will give here a brief summary of how can relate our quantities to the standard Bardeen potentials. A detailed 
analysis of the connection between these formalism is given in . Here we limit ourselves to give the main results 
for scalar perturbations. 

In Bardeen's approach to perturbations of FLRW spacetime, the metric gab is the fundamental object. If gab is the 
background metric and gab = gab + 5gab defines the metric perturbations Sgab in these coordinates. 
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The perturbed metric can be written in the form 

ds^ = S\f]){-{1 + 2A)drf - 2B^dx°'dri + [(1 + 2Hl)1c.p + 2Hlp\dx°' dxf^} , (Bl) 

where rj is the conformal time, and the spatial coordinates are left arbitrary. This spacetime can be foliated in 
3-hypersurfaces E characterized by constant conformal time r/ and metric jaf3- 

The quantities A and are respectively the perturbation in the lapse function (i.e. the ratio of the proper time 
distance and the coordinate time one between two constant time hypersurfaces) and in the shift vector (i.e. the rate 
of deviation of a constant space coordinate line from the normal line to a constant time hypersurface) , represents 
the amplitude of perturbation of a unit spatial volume and is the amplitude of anisotropic distortion of each 
constant time hypersurface (3^ . 

The minimal set of perturbation variables is completed by defining the fluctuations in the energy density: 



/i = /i + (5/i , 6 — 5^i/^, (B2) 

and the fluid velocity: 

u" = u'' + Su" , 5u°' = , 5u^ = -u°A , (B3) 

together with the energy flux and the anisotropic pressure i^ab which are GI by themselves. 

These quantities are treated as 3-fields propagating on the background 3-geometry. With suitable choice of boundary 
conditions [28*1, these quantities can be uniquely (but non-locally) decomposed into scalars, 3- vectors and 3-tensors: 



B^ = B\^+Bi, (B4) 

-S I prTT 

T(a\l3) + -"Ta^ > 



where the slash indicates covariant differentiation with respect to the the metric of E. In this way Vapf — 
/|/3q ~ 5^^/ ^^"^ — f^'^h is Laplacian. The superscript S' on a vector means it is solenoidal {Ba^°' — 0), and 
TT tensors are transverse (ffj^J^i^ = 0) and trace-free. 

On the base of (jB4p and (jB5[) . it is standard to define scalar perturbations as those quantities which are 3-scalars, 
or are derived from a scalar through linear operations involving only the metric 7^/3 and its | derivative. Quantities 
derived from similar operations on solenoidal vectors and on TT tensors are dubbed vector and tensor perturbations. 
Scalar perturbations are relevant to matter clumping, i.e. correspond to density perturbations, while vector and tensor 
perturbations correspond to rotational perturbations and gravitational waves. 

Given the homogeneity and isotropy of the background, we can separate each variable into its time and spatial 
dependence using the method of harmonic decomposition. In the Bardeen approach the standard harmonic decompo- 
sition is performed using the eigenfunctions of the Laplace-Beltrami operator on 3-hypersurfaces of constant curvature 
S (i.e. on the homogeneous spatial sections of FLRW universes). In particular these harmonics are defined by 

y2yjfc) ^ „fc2yjfe) ^ (g7) 

V'Y^f = -k'Y^f , (B8) 

where y«,ri'=\ri^) 

are the scalar, vector and tensor harmonics of order k. In this way one can decompose scalars, 

vectors and tensors as 

/ - f{v)Y (B9) 

B^ = i?(0)(r7)ri")+i?(i)(r,)yii) , (BIO) 

Ht^, = H^"\n)Y^"^ + H^r'\v)Y^'J + K'Hri)Y^f ■ (Bll) 

The key property of linear perturbation theory of FLRW spacetimes, arising from the unicity of the splitting (|B4p . 
(|B5I) . is that in any vector and tensor equation the scalar, vector and tensor parts on each side are separately equal. 

All the quantities define above can be decomposed in this way. However, before proceeding, one should note that the 
quantities A, Ba, H^, H^p, S, v" change their values under a change of correspondence between the perturbed "world" 
and the unperturbed background, i.e., under a gauge transformation. In order to have a gauge-invariant theory one 
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has to look for combinations of these quantities which are gauge invariant. Bardeen constructed the following GI 
variables to treat scalar perturbations ^15] (giving only their fc-space representation); 



Vs 



S' 

s 

kS 



Y 



Bio) 



(0)/ 



Y 



-H, 



(0)/ 



y(0) 



S' 



Y 



(B12) 
(B13) 
(B14) 
(B15) 



where the prime denotes derivative with respect to the conformal time rj. Note that there is not a preferred choice of 
GI density perturbation in this context, as many other GI combinations are possible [33] . 

The variables covariantly defined in the main text are, by themselves, exact quantities (defined in any spacetime) 
and are GI by themselves, therefore, to first order, we can express them as linear combinations of Bardeen's GI 
variables. In [16] this expansions is given in full generality. Here we will limit ourselves to a few examples, giving only 
the scalar contributions and refer the reader to [16] for details. 

The scalar part of the shear, electric part of the Weyl tensor, energy flux and anisotropic pressure are given by 



- iA:2($^ - 



s 



(0) 



S 



Tq/3 



52 



(0) 
a(3 



while the scalar parts of the energy density, expansion and 3-curvature scalar gradients can be written as 

I?„ - -kae„,{ri)Yi°^ , 

^ |-3fc (<i>i, - 1$^) + [UK e) Inhs'] } , 
c„ = -45fc(fc2-3if)(^$^-^yi°))yi"). 

Finally the key covariant scalar perturbation variables are given by 

A^^f^^V^Pa = -fc2e„(77)r, 
C = 5(^V"Ca = k^S^R*{r])Y. 
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The relations above can be used to give an intrinsic physical and geometrical meaning to Bardeen's variables, and also 
to recover his equations. For example, from the (|B16p the variables Vs can be recognized as the scalar contribution 
to the shear. 

The variables Sm (|B15I) . interpreted by Bardeen as the usual density perturbation S^/ ^ within the comoving gauges 
V — B = 0, acquire a covariant significance as the scalar "potential" for the fractional density gradient Va (|B2ip . 
Obviously, this quantity also represents the potential for the divergence A (jB23|) of T>a (or its harmonic component). 

The two independent GI metric potentials and can be combined in such a way to give 
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the former $^ is a stress potential while the latter $jv plays exactly the role of a Newtonian gravitational potential. 
This last interpretation follows directly through (jB17[) . where the scalar part of Eap has exactly the same form it has 
in Newtonian theory Ea/s — ^a/s^N i2fi], independently of any gauge choice. For a perfect fluid {nab — 0) and 
are proportional to each other i.e. = — 1 however in our case, the "curvature" fluid is imperfect, so this is not 
the case. In fact, we can see from equation ([M)) that the anisotropic pressure iTab is related to the shear of the matter 
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flow and therefore the scalar potentials $a and are related to the Bardeen shear variable Vg°^ This relationship 
can be calculated explicitly using ([M)) together with equations (|B16|) . (jB19[) and the expression for the Ricci scalar 

R = 2 (^Q + + ^Rj . After a lengthy calculation we obtain 
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APPENDIX C: COVARIANT IDENTITIES 

On a flat Friedmann background, the following covariant linearized identities hold: 
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where Va = V^a) and W^ofc = W^(a&) are first order quantities. 
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